import gzip
import  Bio
import os,sys
import pdb

from Bio.PDB.PDBParser import PDBParser
from Bio.SeqUtils.ProtParam import *




pdb_stock='/net/stockage/bioinfoDM/pdb/'# repertoire de stockage des fichiers pdb
#path_xml='/net/cremi/nfontain/vrac/xmltest_pdb.xml'# fichier de sortie 


def compute(outputPath):
	path_xml=outputPath        # fichier de sortie 
	list_pdb=sorted(os.listdir(pdb_stock)) # tri de la liste des chemins des pdb

	#list_invader est la liste des molecules qui ne sont pas strictement des proteines reperees par script bash adnCounter.sh
	list_invader = ["4dnb.gz","4bna.gz","3ltu.gz","3k1v.gz","3hg8.gz","3dix.gz","363d.gz","2v6w.gz","2pn3.gz","2b8s.gz","290d.gz","266d.gz","236d.gz","1zx7.gz","1z43.gz","1y0q.gz","1qpy.gz","1nlc.gz","1mme.gz","1d82.gz","1d61.gz","1d26.gz","1ana.gz","159d.gz"]
	list_invader.sort
	list_aa = { "ALA" : "A", "ARG" : "R", "ASN" : "N","ASP" : "D", "CYS" : "C", "GLU" : "E","GLN" : "Q", "GLY" : "G", "HIS" : "H","ILE" : "I", "LEU" : "L", "LYS" : "K", "MET" : "M", "PHE" : "F", "PRO" : "P","PYL" : "O", "SEC" : "U", "SER" : "S","THR" : "T",  "TRP" : "W", "TYR" : "Y","VAL" : "V"  }
	for i in list_invader:
		list_pdb.remove(i)


	#ouverture du fichier de sortie
	output = open(path_xml ,"w")
	output.write("<?xml version='1.0' encoding='UTF-8'?>\n"+"<database>\n")

	parser = PDBParser(PERMISSIVE=1) #(PERMISSIVE=1) gestion d'erreurs lies aux pdb


	# parcours des fichiers pdb
	for k in range (len(list_pdb)): 
		pdbfile=list_pdb[k]
		invader2=True # validite de la molecule(ne contient pas d'acides nucleiques)
		seq = ''
		nb_chain=0
		percent_cys=0
		nb_cys=0
		nb_sheets=0
		nb_helices=0
		nb_turns=0
		nbre_aa = 0
		seqs=[]
		seqtotal=""
		pI=0
		decompress = gzip.open(pdb_stock+pdbfile, 'rb')  #decompression
		flag = 0
		while flag == 0 :
			l=decompress.readline()

			if not l=="":
				k=[]
				k=l.split()
				if k[0]=="SHEET":
					nb_sheets= nb_sheets+1
				if k[0]=="HELIX":
					nb_helices= nb_helices+1
					
				if k[0]=="TURN":
					nb_turns= nb_turns+1
				
			else :
				flag = 1
		decompress.close()
		decompress = gzip.open(pdb_stock+pdbfile, 'rb')
		pdb_id =   pdbfile.split('.')[0]
		structure = parser.get_structure(pdb_id ,decompress)
		header = parser.get_header()
		model = structure[0]	
		name = header['compound']['1']['molecule']
		#parcours de la molecule
		for chain in model :
			nb_chain=nb_chain+1
			seq=''
			#parcours de la chaine
			for residue in chain :
				if (residue.resname in [" DA", " DT", " DG", " DC","  A", "  T", "  G", "  C", "  U","UNK"]):
					invader2 = False 
				else :
					for i in list_aa.keys():
						if residue.resname == i :
							seq= seq +list_aa[i]
							seqtotal= seqtotal +list_aa[i]
							nbre_aa = nbre_aa + 1
			seqs.append(seq)

		if (invader2 == True ) :
			#calcul via la classe Bio.SeqUtils.ProtParam 
			prot=ProteinAnalysis(seqtotal)
			pI=prot.isoelectric_point()	
			aminoDict=prot.get_amino_acids_percent()
			percent_AA_hydro=aminoDict['F']+aminoDict['L']+  aminoDict['I']+ aminoDict['M']+ aminoDict['V']
			percent_cys=aminoDict['C']
			nb_proline=aminoDict['P']*nbre_aa
			
			#generation du xml
			output.write("<molecule>\n")
			output.write("<header>" + name +" "+pdb_id+"</header>\n")          
			output.write("<size>" + str(nbre_aa)+"</size>\n")
			output.write("<nb_chain>" + str(nb_chain)+"</nb_chain>\n")
			output.write("<sequences>\n")
			for s in seqs : 
				output.write("<sequence>" + str(s)+"</sequence>\n")
			output.write("</sequences>\n")
			output.write("<cysteines>" + str(percent_cys)+"</cysteines>\n")
			output.write("<sheets>" + str(nb_sheets)+"</sheets>\n")
			output.write("<helices>" + str(nb_helices)+"</helices>\n")
			output.write("<turns>" + str(nb_proline)+"</turns>\n")
			output.write("<pI>" + str(pI)+"</pI>\n")
			output.write("<percent_of_AA_hydrophobes>" + str(percent_AA_hydro)+"</percent_of_AA_hydrophobes>\n")
			output.write("</molecule>")
	output.write("</database>")

